{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 23\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Code from the previous chapter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "degree"
      ],
      "text/latex": [
       "$\\mathrm{degree}$"
      ],
      "text/plain": [
       "<Unit('degree')>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = UNITS.meter\n",
    "s = UNITS.second\n",
    "kg = UNITS.kilogram\n",
    "degree = UNITS.degree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>x</th>\n",
       "      <td>0 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>y</th>\n",
       "      <td>1 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>g</th>\n",
       "      <td>9.8 meter / second ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mass</th>\n",
       "      <td>0.145 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diameter</th>\n",
       "      <td>0.073 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho</th>\n",
       "      <td>1.2 kilogram / meter ** 3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C_d</th>\n",
       "      <td>0.3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>angle</th>\n",
       "      <td>45 degree</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>velocity</th>\n",
       "      <td>40.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>20 second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>0.2 second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "x                             0 meter\n",
       "y                             1 meter\n",
       "g             9.8 meter / second ** 2\n",
       "mass                   0.145 kilogram\n",
       "diameter                  0.073 meter\n",
       "rho         1.2 kilogram / meter ** 3\n",
       "C_d                               0.3\n",
       "angle                       45 degree\n",
       "velocity          40.0 meter / second\n",
       "t_end                       20 second\n",
       "dt                         0.2 second\n",
       "dtype: object"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_end = 20 * s\n",
    "dt = t_end / 100\n",
    "\n",
    "params = Params(x = 0 * m, \n",
    "                y = 1 * m,\n",
    "                g = 9.8 * m/s**2,\n",
    "                mass = 145e-3 * kg,\n",
    "                diameter = 73e-3 * m,\n",
    "                rho = 1.2 * kg/m**3,\n",
    "                C_d = 0.3,\n",
    "                angle = 45 * degree,\n",
    "                velocity = 40 * m / s,\n",
    "                t_end=t_end,\n",
    "                dt=dt)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params):\n",
    "    \"\"\"Make a system object.\n",
    "    \n",
    "    params: Params object with angle, velocity, x, y,\n",
    "               diameter, duration, g, mass, rho, and C_d\n",
    "               \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    angle, velocity = params.angle, params.velocity\n",
    "    \n",
    "    # convert angle to degrees\n",
    "    theta = np.deg2rad(angle)\n",
    "    \n",
    "    # compute x and y components of velocity\n",
    "    vx, vy = pol2cart(theta, velocity)\n",
    "    \n",
    "    # make the initial state\n",
    "    R = Vector(params.x, params.y)\n",
    "    V = Vector(vx, vy)\n",
    "    init = State(R=R, V=V)\n",
    "    \n",
    "    # compute area from diameter\n",
    "    diameter = params.diameter\n",
    "    area = np.pi * (diameter/2)**2\n",
    "    \n",
    "    return System(params, init=init, area=area)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drag_force(V, system):\n",
    "    \"\"\"Computes drag force in the opposite direction of `V`.\n",
    "    \n",
    "    V: velocity Vector\n",
    "    system: System object with rho, C_d, area\n",
    "    \n",
    "    returns: Vector drag force\n",
    "    \"\"\"\n",
    "    rho, C_d, area = system.rho, system.C_d, system.area\n",
    "    \n",
    "    mag = rho * V.mag**2 * C_d * area / 2\n",
    "    direction = -V.hat()\n",
    "    f_drag = direction * mag\n",
    "    return f_drag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Computes derivatives of the state variables.\n",
    "    \n",
    "    state: State (x, y, x velocity, y velocity)\n",
    "    t: time\n",
    "    system: System object with g, rho, C_d, area, mass\n",
    "    \n",
    "    returns: sequence (vx, vy, ax, ay)\n",
    "    \"\"\"\n",
    "    R, V = state\n",
    "    mass, g = system.mass, system.g\n",
    "    \n",
    "    a_drag = drag_force(V, system) / mass\n",
    "    a_grav = Vector(0, -g)\n",
    "    \n",
    "    A = a_grav + a_drag\n",
    "    \n",
    "    return V, A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func(state, t, system):\n",
    "    \"\"\"Stop when the y coordinate is 0.\n",
    "    \n",
    "    state: State object\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: y coordinate\n",
    "    \"\"\"\n",
    "    R, V = state\n",
    "    return R.y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Optimal launch angle\n",
    "\n",
    "To find the launch angle that maximizes distance from home plate, we need a function that takes launch angle and returns range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def range_func(angle, params):  \n",
    "    \"\"\"Computes range for a given launch angle.\n",
    "    \n",
    "    angle: launch angle in degrees\n",
    "    params: Params object\n",
    "    \n",
    "    returns: distance in meters\n",
    "    \"\"\"\n",
    "    params = Params(params, angle=angle)\n",
    "    system = make_system(params)\n",
    "    results, details = run_ode_solver(system, slope_func, events=event_func)\n",
    "    x_dist = get_last_value(results.R).x\n",
    "    print(angle, x_dist)\n",
    "    return x_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's test `range_func`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "45 102.67621633303261 meter\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "102.67621633303261 meter"
      ],
      "text/latex": [
       "$102.67621633303261\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "102.67621633303261 <Unit('meter')>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "range_func(45, params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And sweep through a range of angles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "20.0 79.83369234082717 meter\n",
      "23.0 86.20478395495141 meter\n",
      "26.0 91.4709452133534 meter\n",
      "29.0 95.83255039259141 meter\n",
      "32.0 99.10671362732461 meter\n",
      "35.0 101.51073514079798 meter\n",
      "38.0 102.9023803583573 meter\n",
      "41.0 103.40314791198459 meter\n",
      "44.0 103.02100653301224 meter\n",
      "47.0 101.73693469942923 meter\n",
      "50.0 99.56858789440345 meter\n",
      "53.0 96.55360348120037 meter\n",
      "56.0 92.7050051463521 meter\n",
      "59.0 88.04022482508677 meter\n",
      "62.0 82.58080180144869 meter\n",
      "65.0 76.32548350395584 meter\n",
      "68.0 69.32998156572299 meter\n",
      "71.0 61.63270335890529 meter\n",
      "74.0 53.24654950607608 meter\n",
      "77.0 44.252471292784826 meter\n",
      "80.0 34.68860100116735 meter\n"
     ]
    }
   ],
   "source": [
    "angles = linspace(20, 80, 21)\n",
    "sweep = SweepSeries()\n",
    "\n",
    "for angle in angles:\n",
    "    x_dist = range_func(angle, params)\n",
    "    sweep[angle] = x_dist"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the `Sweep` object, it looks like the peak is between 40 and 45 degrees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap23-fig01.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3ddXhU19bA4d9ESUhw16ALLe4UaKm7G6XUjcqt3dtb70fd21ul7k69lJYKUopD8YW7NjghxOb7Y5+0Q5rAQJOcmWS9z5MnmaPrzEzOOnufffYOBINBjDHGmEgT43cAxhhjTGEsQRljjIlIlqCMMcZEJEtQxhhjIpIlKGOMMRHJEpQxxpiIFOd3AKb0iMhyoHGByTuAmcDNqjq5tGOKdCISA7wGnAVsUtWC719x7+9UYKqqrhKRAcDPQKqq7izJ/RYRy6HAG0B94F+q+mKB+UHgRFX9urRjKxDHPcAJqtrVzzjyRVo80cxKUOXPbUBd76ceMBDIAr4VkRQ/A4tQvYAhwOlAn5LckYg0BkYAlb1JE3Cf066S3O8+3AUsBFoB7/kUgynHrARV/uxQ1fUhr9eJyIXAKuBw4EtfoopcVbzf36lqST/VHgh9oapZwPoili0NVYBvVXW5jzGYcswSlAHY4/3OBRCRisCjwClATdxJ8iVVvc+b/waQASQBZwDpwCsh8wPAvcAVQAVcFVlH4E1VfcNb5mbgWqA6f1UxTiwsuDDiaQs8C3QDdgNfANer6t9KHiIS58U2CFd1tRn4ALhRVXMLLHsh8Lr3Mk9E7vX+3qv6RkR+wVXL3eytcw3wEXADEA98B1yRH4+InAHcCbQElgK3qeoXwDJvk7O9ff1CSBWfiNQFHgGO9t7X77zjXOdtNwhc5L2vrYHZuKq534p4X4vcXkh1cFcRuUtVA4VtI2Rb+/uM/nyPQtb5s4owjO9ULK70fxlQDZgCXKOqc73NxYrII8AluPPaB8C1XpIvLN4bgKuANGAn8A1wlfc+X8j+P8OzgGFAI+BH3OdYSVUvLGRfLYFngH7AJuBD4E5V3VNwWbM3q+Ir50SkOvAy7oQyzpv8BK5q6xRAcP9cw0SkS8iqlwKrga7AKwXm/wf3z38prlpMgP4h+7wCuB64GugEfAv8JCJNighzf/G8ByhwCHAirtry1iK2dQswGFdt18J7PdTbdkEfesuCq2p7rIhtFnQIcChwBO6EehpwJYCIHO5t922gPTAc+EhE2gDdvfUHFNyXiMTjToSNgONwpd36wOfeBUG++3DJrxeu6nZ4YQGGsb1uwO/A496x708435n92dd36i7gOuBfuO/MKuAbL3GBuwCqBvTEfWYXeT9/IyLnAvcAN+K+AxcCJ+MuqPLt6zPsDbwLPO/tdxYuoRW2rwrAKFwC6+TFdgzwdBjvR7lnJajy53ERecj7OwaIxSWmI1R1uzf9V9zV73Tv9aMichfQBpjmTVuiqnd6f9/nXZF28eZfCwxT1a8AROR83Akl323Arar6jff6Aa9BwFDgZv5uf/Gk4aomV6jqUhE5CXdyLsxc4EJVHeO9Xi4it3jb+jR0QVXdLSJbvb/Xe8dSxGb3Eg9c5q0zV0S+w7034BL3l6qan4Ce9u79JeOurgHSvSv50G0eDTTHfU5rvVjOxpW6jgB+8JZ7VlW/9eY/CnwhIomFXK3vc3uq+oOI5AA7C1QJFyWc78z+FPqdEpHpuPftPlUd4cU6FJdkqnrLbwGuVNUcYJGIjAU6FLGftbjvQH7jjhUiMsaLNd++PsNrgK9UNT/J3OZdeBTmXCAbGOpVEauIXAmME5F/h/zPmUJYgip/HgTeARJxJZjTgHtCqkrw5p8gIoNx1VAdgRRcMsu3qMB2dwDxIlID1/hiSv4MVd0sIgsBvJNxI+BlEXkpZP1E/qpqLGh/8fwHeA64SkRGAR+r6ueFbUhVvxSRw7zqoJa4K+W0Asf2TxW8z7cdqOj93QZ39R0a0/0AIpK2j222xSXgtSHrrfaq4tryV4JaWGC/4P7PC7634W4vXOF8Z/an0O8UUANXbRj6ndoB3AR/XjSs8JJTvq246sK/UdUxItJVRO7DNQBpiyv1vRW67318hodQ4DMEfuOvZBmqLdAU2BFywRHAXRy2IPzkXS5ZFV/584eqLlbVuao6FFd6+kpEGoYs8xqu+mI37p+2F+4fPlRhJZQA7moRiv5u5Z+whuBOYvk/rfGqUAqxz3i85s+NcfeWquCqzF4ubENeE+ARXhyfAsfiqrLCVVhDiYIXekW9N/nzDqaxxe4ipuef7MLZ98FsL1z7+87sdczevcCCioo9f/q+3rfcQqYVet/Mu8c0DlclOBI4h783DtrX+5hN+O9RHC55hX7XO+CS07wwt1FuWYIyVwM5wAsAIpKKqye/QFVvU9UPcf+slSniHz6Uqm7D3Uf4896DiFTG/UPmz18P1PcS5WJVXYyrNjm64Pb2F4+IVBKRZ4Ggqv5PVY/3tjWoiBCHAjep6k2q+jawHJfc9ntsnvx958cXAIq6d1aYhUDn0AkiMkpE/sW+T8DzgcZew4b89erhYl9wAPsv9u2F+Z3Z633DlSrC4n1nNhLyvolIkohsEJEeBxKrZyjwqKperaqv4u4htSD878AcQr7fnm5FLDvf2/bqkO96NeAhIOGAIy9nrIqvnFPVP0TkVuAVETkZ12BhF3CaiCzDVdc9ivvnTQxzs08Cd4jIUtzN4ftw1SP5J+BHgLtEZB2u2uY8XFLpX8i2MvcVj6pu9+5f1RORO7x1TiGkOqiAdOB47x5FJf66jxHusU3B3TO7Cvged7+tWpjrAjyFu/9wDe7q/Xhc665/4VqTAXQUkVUF1huNK+l9ICI3etOewCW80Qew/5LY3j4/I2+ZKcD1IvIWsM3b14G0YnsSuNOrglwI3O5tZwauFHwg0oHDvIYpAVxLvTa4Vo/heBr4TUSuw7XuOxvXGKhgFSW4qs+7gDdFZBiQimsAstJLvGYfrARlwFXP/Ir7x4vHJYwjcFUQb+FOxN/y96vGojyF64HgTVz1xmJgBX9VmzyNa6X2iLePs4AzVPXXghtS1eww4jkN18hgAjAJd8I8r4jYLsSVeGYDn+NKUK8ewLH9iLv6vQ93/yAP16Q5LF6T7wtwCXkucDFwiqrOV9V0L5ZXcNWVoesFcYl3E675+Y+4m/0Di2pKvZ84im17YX5Gj3vTRgFf407cqw9gN4/y13szA9ey8PiDOXZcC9IgMBWXjBNx92Y772ulfKo6FddC8AZc6asr7tGGv8XiNUs/CncRNBn4ir8uysx+BGxEXVPcROQ4YHpIy7c44A/gJFUd62twxvxDXrXiztCGRSLyDTBJVf/Pv8jKHqviMyXhEuBGr5lwJq76aiuudGNMtOsO/FtEBuFqBo7GPUf2b1+jKoMsQZmScA2uZ4exuO/Yr8BR9uS8KSOexz2a8AGu6m4+rop67r5WMgfOqviMMcZEpDJfghKRRFwT0HUU/qyEMcYY/8TiGr1MKVjLUuYTFC45jdvvUsYYY/x0KDA+dEJ5SFDrAN59913q1KnjdyzGGGNCrF+/nkGDBoF3rg5VHhJULkCdOnVo0KCB37EYY4wp3N9uwdiDusYYYyKSJShjjDERyRKUMcaYiGQJyhhjTESyBGWMMSYiWYIy5UZObg5Zudn7X9AYExHKQzNzU0bl5uWyI2sX2zN3sH3PDrbt2cH2zJ3eb+/1np1//p2R7QaRjY+JIzkhmZT4ZComJFMxIYmK8ckkJySRkpBMcnwyKQnevPgkKiZU9H4nkxRfgZiAXdcZUxosQZmIlpuXy4qtq1nwxxIWpS9jy+5tbN/jktDOPbsIHsDo6TGBGAKBANl5OWzL3M62zO0HHE9sTCxpVRrQonoTWlRrQosaTahdsQaBQLiDsRpjwmUJykSUzJw9LE5fxoI/lrBg0xIWpi8lM6fwTtADBEhNTKFSYgqVE1OpVCH1z78rV0ilUqL7qexNr5iQTIAAWbnZ7MrKYFd2hvd7t/vtTduZlUFG1m52ZmeQkbX3MrtzMlmyeQVLNq/gO34BIDUxhRbV0lzSqt6E5tXSSE5IKsV3zZiyyRKU8dX2zB1eMlrMgj+WsGzLSnKDeXstUyelJq1qNEdqNKV2Sg2XeCqkkppQkdiY2APeZ2JcAolxCVSjygGvm5G9myWbV7AofdmfP9v37GT6ujlMXzcHcImzfqU6XsJyiathpXrExFjVoDEHwhKUKTXBYJANu/5wychLSGt3bNhrmUAgQNOqjWhVoxmtajanVY1mVEmq7FPEf5ccn0T72q1oX7sV4I5p464/WJS+nEXpy1icvoylW1exevs6Vm9fx8/LJgBQIS6RZtUa06J6EzrWaUOrms3tXpYx+2EJypS45VtW8f2ScUxbM4stmdv2mpcQG0+L6k1oXbM5rWo0p0X1JiTFV/Ap0gMXCASonVKT2ik16du4GwBZudks37KKxZuXs9ArZW3alc7cjQuZu3Ehn88fRc3kahya1oN+aT2ol1rb56MwJjJZgjIlYk9OFr+tmsYPi8eyaPPyP6enJqa40lGN5rSu2Zy0qg2JO4hqukiWEBtPyxpNaVmjKcd507Zmbv/z3tqEldPYlLGZEfNGMmLeSFpUb0L/tB70btiVlMSKvsZuTCTxJUGJSHfga1Wt5b1OwA0RfgauR9snVPXBkOXPAh7ADWo1BrhQVTeWeuBmv9ZsX88PS8YxZtlv7PKadSfHJ9E/rSeHNelN4yr1y2WLtyoVKtG1fge61u/AeYecwvxNixmzfCITV03/817W6zM+pku99vRP60mnOm2Ji7XrR1O+lep/gIgEgEuAxwrMuhcQoBlQGfhORNao6lsi0gZ4FTgWmAo8DHwAHF5qgZt9ysnNYfKamfywZBxzNy78c3qzao05slk/+jTqSmJcgo8RRpaYQAxta7Wkba2WXNL5HKasmcmY5ZOYtWE+k1fPZPLqmaQmVKRPo270S+tBs2qNy2VSN6a0L9HuBY4H7gPuCJk+BFcq2gJsEZHHgCuAt4Dzga9UdTyAiPzXW6aFqi4q1ejNXjbu/IPRS8fz89IJbNuzA4DE2AT6Nu7Okc360rRaY58jjHyJce796tu4O5t3b2X8iimMWT6RVdvW8t3iX/hu8S/UT61Dv7QeHJrWnRrJ1fwO2ZhSU9oJ6kVVvUtEBuRPEJEquKq7eSHLLQDae3+3wZWcAFDVDBFZ5c23BFXK8vLymL5uNj8sGcfMdfP+fFC2YeV6HNnsUPo17mHPAB2kaklVOKnVkZwoR7Bi62rGLJ/E+BWTWbNjPe/P/oIPZn9J21otGdCkF70bdS1z9+6MKahUE5Sqri1kcor3OyNkWgaQHDI/g72FzjelYPPurfy0dAI/Lh1PesYWAOJi4ujVsDNHNuuH1Ghq1VDFJBAIkFa1IWlVGzKow6nMWj+fscsnMmXN78zZqMzZqHw89xvOaHMcfRt3O6hnwYyJBpFwF3aX9zv0sjsZ2Bkyv+Aleeh8U4J2Zu3is3nfMXLRL+Tk5QDuwdkjmh3KgCa9qJSYsp8tmH8iLiaWzvXa0bleO3ZlZfDbqul8pT+wbsdGnpv8Jp/N/44z2h5P70Zd7LkqU+b4nqBUdYuIrMc1kljjTW7FX1V+87x5AIhIMtCIvasETTHLyc3h+yVj+WTut+zMctcQ3et35Kjm/WhXW+xk6IOKCckc0awvhzXpxbgVk/l07res3bGBZya+xmfzRnJmuxPo3qCjfTamzPA9QXneBu4WkVm4Kr2bgae9ee8B4737Vr8BDwIzVHVhYRsy/0wwGGTS6hm8O+tzNuzcBEDbWi0Z3OE0a/QQIWJjYhnQpBd9G3dnzLLf+HTeSFZtX8cTE16mcZUGnN3uBLrUO8SqXE3Ui5QEdRfwODAXN0bVcOBFAFWdLSIXe6/rA5OAM32Ks0xb+MdS3p75KZq+FID6qXUY1OFUutRrbye7CBQXE8vAZn3pl9aDn5dNYMS871ixdTWPjH+RZlUbc1b7E+hYp619diZqBYLB8IcriEYikgYs+/HHH2nQoIHf4USkDTs38d6sL/ht1TQAKiWmcFa7ExjYtK/dgI8iWbnZjF4yjs/mj/pzKJGW1ZtydvsTaVdLLFGZiLR69WoGDhwI0ERVl4fOi5QSlPHBzj27GDFvJCMX/0JuXi7xsfGc0HIgJ7c+iuR4ayoebRJi4zmu5eEMbNqXUYvH8MWC71mYvpRhvzxN65otOLvdibSp1cLvMI0JmyWocig7N5tRi8fw6byR7MrKIECA/mk9Obv9ifYgaBmQGJfASa2O5MhmhzJy0c98paOZv2kR9/z8BO1rt+LsdifSskZTv8M0Zr8sQZUjwWCQ31ZN571Zn7FxVzoA7WoJgzueTpOqDX2OzhS3pPgKnNbmWI5pPoBvFv7I1wt/ZPaGBczesIDejbpycaezqFQh1e8wjSmSJahyQv9YwlszP2VR+jIAGlSqy/kdTqNTXbuJXtYlJyRxZrsTOLbFYXy9cDTf6E9MWDmV2RsWcFGns+jTqKt9B0xEsgRVxu3JyeLNGR8zeul4ACpXqMTZ7U7gsCa9rQFEOZOSWJFz2p/M4U368NLUd5i9QXlm4muMXzmFy7ucR7XkAx9h2JiSZAmqDFu9bR1PTniZVdvXER8Tx0mtjuKkVkdG1YCApvjVSqnBHf2v5+dlE3hr5qdMXzubGzbdy+AOpzOwaR8rTZmIYQmqDAoGg/y8bAKvTf+QrNxs6qXW5l+9LiWtqjWzN04gEODwpn3oWLctr0x9n6lrZzF86rtMWDmVK7oNonZKTb9DNAbrE6WMycjezTMTX+PFKe+QlZtN/7SePHTkrZacTKGqJVXhlr5X8q9el5CamMKcjcrN393HN/ojeXl5fodnyjkrQZUhSzev4MnfXmXDzk0kxiVyWZdz6ZfWw++wTIQLBAL0btSVdrWEN2Z8zPiVU3hz5idMWDWNq7oNpkHlun6HaMopS1BlQDAYZOSin3n79xHk5uXSuEoDbuh1CfUq1fE7NBNFKlVI5bpeF9O7UVdemfY+i9KX8e/vH+D0NsdycuujbfwpU+osQUW5nXt28fzkt5i6dhYARzfvz+COp5MQG+9zZCZada1/CG1qtuDt30fw49LxfDjnKyaunsFV3QbTtFojv8Mz5YglqCi2YNNinp74GukZW6gYn8SV3QfTo0Env8MyZUByQhJXdBtEn0ZdeGnKu6zYuprbRj/MSa2O5Iy2x9sFkCkV1kgiCuXl5TFi3kju+flJ0jO20KJ6Ex4++nZLTqbYtavdikePuYPjWw4kGAzy+fxR3DLqPhZsWuJ3aKYcsBJUlNm6exv/m/QGszcsAODkVkdxdvuT7P6AKTEV4hIZ0ukMejXszAtT3mbN9vXc/fPjnH/IaZwgA+25KVNirAQVRWatn88to+5n9oYFVEpM4bZ+1zKow6mWnEypaFmjKY8cdRsntTqKYDDI279/ylO/vUpmdqbfoZkyykpQUSAnL5eP5nzFF/O/J0iQdrWEa3teRNWkyn6HZsqZ+Nh4zu9wKi2rN+G5SW/y26pprN6+jpv7XEHd1Fp+h2fKGEtQEW7Hnp089utLzN+0mEAgwNltT+TU1scQE2OFX+Of7g06Ur9SHR4b/xKrtq3lvz88xLU9L6JLvfZ+h2bKEDvLRbBNu9K568fHmb9pMVWTKnP3gBs4ve1xlpxMRKhfqQ73H/lvutfvSEb2bh4e9zwfzfmavKD1QGGKh53pItTyLau4Y/SjrNmxnkaV6/PgEbfaaKgm4iTHJ3FTn8s5t/3JBAjwydxveGTcC+zKyvA7NFMGWIKKQLM3LODun55gS+Y22tZqyf8dfpMNhWAiViAQ4NQ2x3Bb/2tISajI9HVzuPWHh1i5dY3foZkoZwkqwoxfMZkHxj7L7pxMejfswm39riE5IcnvsIzZrw512riOias0YMPOTdw++hF+XTnF77BMFLMEFSGCwSBfLvieZya+Tm5eLse3HMh1vS4m3p7YN1GkVkoNhg28hX6Ne7AnN4unf3uNt2Z+Sm5ert+hmShkCSoC5AXzeHPGx7zz+2cAXNDxdIZ0OoOYgH08JvokxiUwtMcQLu58NrGBGL7W0dw35hm2ZW73OzQTZewM6LOs3Gye+u1Vvl30M7ExsVzf62JOkCP8DsuYfyQQCHBMiwHcfdiNVKlQibkbF3Lr9w+xOH2536GZKGIJyke7sjJ4YMz/mLhqOknxFbi93zX0adTN77CMKTatajbjoaP+i1RvSvruLdz10+P8tPRXv8MyUcISlE/SM7Zw14+PMW/TIqpWqMy9h91Eu9qt/A7LmGJXLakKdx92A0c3709OXg4vTnmH4VPeJTs32+/QTISzBOWDlVvXcMfoR1m1fR31K9XhviNusSHZTZkWFxvHJV3O4eruFxAfG8/opeN5eNwL1o+f2SdLUKVs3sZF3P3T46Tv3oLUaMaww2+mZsXqfodlTKkY0KQXww6/icqJqczaMJ9hvzzNjj07/Q7LRChLUKVo4qrp3DfmGXZl76Z7/Y7c2f86UhIr+h2WMaWqabXG/N/Am6mZXI1Fm5dz909PsDljq99hmQhkCaqUjFz4M09OeIWcvByObt6fG3tfRkJcgt9hGeOLuqm1GDbwFhpWqsvq7eu488dHWbdjo99hmQhjCaqE5QXzeOf3z3h9xkcECXLeIadwceezrcNXU+5VS67CvYffRItqaWzK2MxdPz7Gsi2r/A7LRBA7S5agYDDIWzM/5csF3xMbiGFo9yGc0vpoG4HUGE9KYkXuHHA9h9RuzbY9O7jn5yeYt3GR32GZCGEJqgSNmDeSbxf+RGxMLLf0vYr+TXr6HZIxEadCfAX+c+hV9GzYmd3Zmdw/9n9MXTPL77BMBLAEVUJGLRrDh3O+IhAIcH3Pi+lcr53fIRkTseJj4/lXz0s4omlfsnOzeezXlxi7fJLfYRmfWYIqAeNXTOG16R8CcFmX8+jZsLPPERkT+WJiYris63mc0vpo8oJ5PDvpDb5d+JPfYRkfWYIqZtPXzuG5SW/82SDiiGZ9/Q7JmKgRCAQ475BTGNzhdADemPExH87+imAw6HNkxg+WoIrRgk2LeWLCcHKDeZzU6ihOaX203yEZE5VObHUEV3UbTCAQ4NN53/Lq9A9sKPlyyBJUMVm+ZTUPjXuerNxsDm/ah0GHnOJ3SMZEtcOa9uam3pcTHxPH94vH8szE18nJzfE7LFOKLEEVg/U7NnL/2P+Rkb2bHg06cXmX86wpuTHFoHuDjvy33zUkxVVgwsqpPDL+BfbkZPkdlikllqD+oc0ZWxnmDcbWvnYrrut5kT2Ea0wxaldbuPuwf5GamMLM9fO475en2Zm1y++wTCmwM+k/sGPPTu4b8wybdqXToloat/S5woZoN6YENK3WmGGH30T15Kpo+lLu+elJtuze5ndYpoRZgjpImdmZPDT2OVZvX0eDSnW5td9QKsRX8DssY8qsepXqMGzgzdRPrcPKbWu456cn2GrDyJdplqAOQnZuNo/++hKLNi+nZsXq3NH/OlITU/wOy5gyr0ZyNe4deBONqzRg3c6N3P/LM1bdV4ZZgjpAeXl5PDPxdWZvWEDlxFTu6H8d1ZKr+B2WMeVGpcQU7uh/LfVSa7Ni2xoeHPucDXxYRlmCOgDBYJDhU99l0uoZJMcncXv/a6mbWsvvsIwpdypXqMQdA66jRnI1FqUv45HxL5JlQ8iXORGToESkp4hMFpFtIrJYRC71pieIyHAR2Swim0Tkv37F+O6sz/lp2QQSYuO59dCrSava0K9QjCn3aiRX464B11OlQiXmbFSenPAyOXm5fodlilFEJCgRiQG+AJ5R1crAucCzItIBuBcQoBnQDRgiIheUdoyfzx/157AZN/a+nFY1m5d2CMaYAuqk1uKO/teRklCRaWtn8/ykN8nLsx4nyoq4cBYSkVigC9AVqAXkAuuBKao6sxjiqOptNyAiASAI5ABZwBDgQlXdAmwRkceAK4C3imG/YRm9ZDzvzfqcAAGG9hhiPZMbE0EaVanPbf2u4f9+eYrxK6dQIS6Ry7raw/JlwT4TlIhUBa4DrgKqA0uBdCAWqAE0FpF1wIvAc6q69WCCUNV0EXkWeBN43dv+9cA6oC4wL2TxBUD7g9nPwZi65ndenvYeABd1Pou+jbuX1q6NMWFqXj2N/xx6NQ+MfZbRS8eTFF+B8zucZkkqyhVZxedVo80AGgGXACmqKqraW1V7qGozoBpwJdAWmCMiQw4mCK+KLxM4D0gCBgB3Ayd5i2SELJ4BJB/Mfg7GmOWTCAaDnNXuRI5pMaC0dmuMOUBta7Xkpt6XExuI4SsdzYh5I/0OyfxD+ypBdQO6qeqmohZQ1e3AN8A3IlIXuB1XCjpQpwF9VPUW7/UYEXkVV70HLmnlSwZ2HsQ+DspFnc/i6Ob9aVurZWnt0hhzkDrXa8e1PS/m6Ymv8uGcr0iKr8BxLQ/3OyxzkIpMUKp67YFsSFXXAdccZBwNgcQC03KATbh7XQKs8aa3Yu8qvxJVLakK1ZLsOSdjokXvRl3IzNnDi1Pe5o0ZH1MhrgKHN+3td1jmIITVSAJARJoCrfl7IkFVR/zDOL4HHhSRy4GXgc7AZcClwErgbhGZBaQANwNP/8P9GWPKsMOb9mZ39m7enPkJL019h6T4RHo17OJ3WOYAhduK7xbgIVzruoJPwwX5h/eEVHWuiJwGDAMexZWablXVL0RkFPA4MBd3z2w4rlGGMcYU6XgZyO6cPXw05yuemfg6ibGJ1gI3yoRbgroFd3/pEVUtkYcMVPVb4NtCpmcCQ70fY4wJ2+ltjiUjezdf62genzCc2/tdQxu7nxw1wn1QNwCMKKnkZIwxJSEQCDC4w2kMbNqX7NxsHh73AovTl/sdlglTuAnqOeBWEUkoyWCMMaa4BQIBLutyLn0adWV3TiYPjH2WlVvX7H9F47twq/g+BsYC54nIemCvkpSqNi3uwIwxprjExMQwtMeFZObsYdra2dw35hn+7/CbqGOdPUe0cBPUO8By4AP2fmjWGGOiQlxMLDf0voyHxj7HnI3KsF+e5r4j/k3VpMp+h2aKEDGjnPUAAB3GSURBVG6CagV0UNWFJRmMMcaUpITYeP7d90qG/fI0izYv55HxL3DvYTeSEGd3LyJRuPegpuAeljXGmKhWIb4C/zn0ampVrM6SzSt4dvKb5AWt/VckCrcE9S7wmoi8DyyhwLNQqvp8cQdmjDElpVKFVP5z6NXc8eOjTFw1nY9Sa3NO+5P2v6IpVeEmqFtx/d+dWMi8IGAJyhgTVRpWrscNvS7jwXHPMmLeSOql1qZfWg+/wzIhwkpQqtqkpAMxxpjS1rFuGy7qdBavTf+QF6e8Q62KNWhVs5nfYRnPvobbOKChM0QkICIX//OQjDGm9BzTYgDHNB9ATl4Oj/76Iht2FjmAgyll+2ok0UlEZonI1d5QGoUSkVoi8i9cX3kdiz1CY4wpYUM6nUGHOm3YsWcnD497gYys3X6HZNj3cBv/EpGuwF3AUyIyD5eE/sB1fVQT6AC0BEbihmWfXPIhG2NM8YqNieWGXpdyx4+Psnr7Op787RVuPfRqYmNi/Q6tXNtnM3NVnaqqJ+GamA8H9gDNgSa4RhPPAE1V9URLTsaYaJackMSth15NpcQUfl8/jzdmfOx3SOVeuI0klmEt9YwxZVytlBrc0vdK7v35KUYtHkP9SnU4psUAv8Mqt8J9UNcYY8oFqdGMq7oNBuD1GR8xc91cnyMqvyxBGWNMAYemdef0NscRDAZ5csIr1vu5TyxBGWNMIc5sdzy9GnZhd04mD49/gW2Z2/0OqdyxBGWMMYWICcQwtPsFNK+WxqZd6Tw2/iWycrP3v6IpNmEnKBFpLCKPiMjnIlJXRC4UkV4lGZwxxvgpIS6Bf/e9kurJVdH0pbw45R2CwaDfYZUbYSUoEemBewaqA3AskIR7KHeMiBTWP58xxpQJVZIq85++V5MYl8j4FZMZMW+k3yGVG+GWoB4F7lPVo4EscA/yAv8HDCuh2IwxJiKkVW3A9T0vJkCAD+d8xYSV0/wOqVwIN0F1wg37XtC7uJ4kjDGmTOta/xAGdzwdgOcmv8mi9GU+R1T2hZug0oEWhUzvBmwovnCMMSZyHd/ycI5o2pfs3GweGf8if+za7HdIZVq4CepZ4CURORfXD19HEbkO17vESyUVnDHGRJJAIMDFXc6hXS1hW+Z2Hvv1JbKtZV+JCStBqepjwAPAg0Ay8AluEMP7gIdLLDpjjIkwcTGx3NjnMmpWrM7SLSt5c8YnfodUZoXdzFxVX1LVNCAVqKKq9VT1KVW1NpfGmHIlJaEiN/W+jLiYOL5fMpZxy62v7JIQVmexInJBIdPADfeeBawBJqmqlXWNMeVC02qNuajTWbw87T2GT32XtKoNaFi5nt9hlSlhJSjgQqAfkAks9Ka1wFX3rQCqAltE5ChVXVzcQRpjTCQ6ollf9I8ljF0xicd/Hc6DR95KUnwFv8MqM8Kt4psJjAIaqmpnVe0MNAC+AD4CauAGLXymRKI0xpgIFAgEuLTruTSsVJe1OzbwkvU0UazCTVAXAbeo6pb8Caq6DbgDuEJVc4GngD7FH6IxxkSuCnGJ3NTncirEJTJh1TRGLR7jd0hlRrgJag+QVsj0JkCu93cikFMMMRljTFSpV6kOV3pjSL058xN7iLeYhHsPajjwuoj8HzAVl9i6AHcCr4hILeBJ4KcSidIYYyJc70Zd0D+WMHLRzzwx4WUePuo2KiWm+B1WVAv3Oai7cAnov8BvwK/Av4GHvGldgG3ANSUTpjHGRL7BHU6jRfUmpGds4X8TXycvmOd3SFEt3BIUqvoQ8JCIVAeyVTV09K6R3o8xxpRbcbFx3ND7Uv7z/YP8vn4eI+Z9xxltj/M7rKgVdoISkU64klI8EPCegwJAVZ8v/tCMMSb61EiuxnU9L+KBMc/y8ZyvaVm9CYfUae13WFEp3PGgbgem4ar0/g3cEvJzc4lFZ4wxUahDnTac0fY4ggR5euJrpGds2f9K5m/CLUFdCtypqveXZDDGGFNWnN7mOBamL+X39fN5csIr3HPYDcTFhl1pZQi/mXl13AO5xhhjwhATE8O1PS+melJVFqYv5Z1Zn/kdUtQJN0F9AgwqyUCMMaasqZSYwg29LyU2JpZvF/7ExFXT/Q4pqoRb3twN/FdEzgQW4Q37nk9VzyruwIwxpixoWaMpF3Q4nddnfMQLk9+mUZX61Eut7XdYUSHcElQy8B4wGdgC7CrwY4wxpgjHtBhAr4Zd2J2TyeO/DmdPTtb+VzLhlaBU9aKSDsQYY8qqQCDAld3OZ8XW1azatpaXp73H0O5DCAQCfocW0Q70Oag2QKw3KYDrf6+Lql5RArEZY0yZkRRfgZv6XM5tPzzM2OWTaFWjOUc06+t3WBEt3AELbweGATuBirhujSp7s78tmdCMMaZsaVi5Hpd1PY9nJ73B69M/pGnVRjSt1sjvsCJWuPegrsANt1EJWAccAtQHJgJTSig2Y4wpc/ql9eDIZoeSnZfDkxNeJiN7t98hRaxwE1Qd4FPv75lAL1Vdj+tVYnBJBGaMMWXVkE5n0qRKQzbs+oM3pn/sdzgRK9wEtQn3sC64Id87eH+vAeoVRyAiUldEPheRbSKyQUSGedMTRGS4iGwWkU0i8t/i2J8xxvglITae63pdTEJsPL8s/82ejypCuAnqC2C4iHQEfgYuEJH+wI3AimKK5Qtc9WFtoCcwRETOA+4FBGgGdPOmX1BM+zTGGF/Ur1SHwR1OB+Clqe+yOWOrzxFFnnAT1M24e03tgK9xAxP+CAzBdRj7j4hID6ApcJ2qZqrqMmAALhkOAe5X1S2quhx4DHdPzBhjotpRzfvRqW47dmVl8NzkN238qALCHbAwQ1WvVNV3VDWoqhcCVXDVfj8WQxxdgNnAPSKyRkSWAKfierCoC8wLWXYB0L4Y9mmMMb4KBAJc1X0wlRJTmL1hASMX/ux3SBFlnwlKRJJF5AQROVZE9hq7WFV3Aseyd/I4WNWAQ4FsXEnqNFyp7SRvfkbIshm4ni2MMSbqValQiSu7nQ/Ae7M+Z+XWNT5HFDmKTFAi0g1YDnwJfAOoiLT25tUXkS9x942K493cA2xX1XtUdY+q/g68gqveA0gKWTYZ9zyWMcaUCV3rd+CIpn3JzsvhmYmvk5Wb7XdIEWFfJajHgN+BhriGC1OAp0WkN646ridwqaoeWgxxLACSRSQhZFocrt+/9bhGEvlaUTylNmOMiRgXdDqDuim1WLltDR/M+sLvcCLCvhJUR+BWVV2jqpuAS4B+wMfAD0BrVX2tmOL4AdeU/XGvWXl7b3/vA28Dd4tIDRFJw1X9vV1M+zXGmIhQIS6Ra3teREwghq8X/sis9fP9Dsl3+0pQqcDq/Beqmg7kAiNU9WzvdbFQ1UygP+7+0zrgO+ARVf0UuAuYA8zFleI+BV4srn0bY0ykaF49jTPbHg/Ac5PfZOee8j1YxP764gsWeJ0HPF8SgajqUuD4QqZnAkO9H2OMKdNOaX00M9fNRdOXMnzqe9zQ+9Jy2+t5uM9BhdpT7FEYY4wBIDYmlmt7XkRSXAUmrp7O2OWT/A7JN/srQV0oIqEt5uKA80Xkj9CFVLVESlXGGFMe1UqpwUWdz+L5yW/x2vQPaV2zObVSavgdVqnbV4JaCVxVYNp6oODghUFKqNrPGGPKq/5pPZm+dg4TV0/nf5Pe4J7DbiA2Jnb/K5YhRSYoVU0rxTiMMcaECAQCXNb1XDR9CfrHEr5Y8D2ntTnW77BK1cHcgzLGGFMKUhNTGNrd9Vfw8ZyvWZy+3N+ASpklKGOMiWCH1GnN8S0HkhvM438TXyczp/y0U7MEZYwxEe7cQ06mUeX6rNu5kbdmfrr/FcoIS1DGGBPhEmLjua7nRcTFxDF6yTimrvnd75BKhSUoY4yJAo2q1Oe8Q04B4IUp77B19zafIyp5lqCMMSZKHNfyMNrXFnbs2ckLU94hGCzY2U/ZYgnKGGOiREwghqu7D6FiQjIz1s3hhyVj/Q6pRFmCMsaYKFI9uSpXdB0EwFszP2XN9vU+R1RyLEEZY0yU6dmwM/3TepKVm80Lk98mLy/P75BKhCUoY4yJQhd2OpOqSZVZmL6UkYt+9jucEmEJyhhjolDFhGQu96r63p/9Bet3bvI5ouJnCcoYY6JUl3rt6du4O1m52bw05R3ygmWrqs8SlDHGRLGLOp1J5cRU5m5cyOgl4/0Op1hZgjLGmCiWmpjCJV3OAeCd30fwx67NPkdUfCxBGWNMlOvZsDM9GnQiM2cPL019t8w8wGsJyhhjyoBLOp9NSkJFfl8/jzHLJ/odTrGwBGWMMWVAlaTKXNjpTADenPExm3dv9Tmif84SlDHGlBGHNu5O57rt2JW9m1emvh/1VX2WoIwxpoxww8SfR1J8BaauncWEVVP9DukfsQRljDFlSPXkqlzQ4XQAXpv+Edszd/gc0cGzBGWMMWXM4U37/Dksx2vTP/Q7nINmCcoYY8qYQCDAFV3PJzEukQmrpjF59Uy/QzoolqCMMaYMqpVSg0HeCLyvTHufnXt2+RzRgbMEZYwxZdRRzfvRqkYztmZu582Zn/gdzgGzBGWMMWVUTCCGK7sPJj42njHLJzJj3Ry/QzoglqCMMaYMq5dam7PbnQjA8CnvkZG12+eIwmcJyhhjyrgTWg6kebU00ndv4Z3fR/gdTtgsQRljTBkXExPDVd0HExsTy+il45m9YYHfIYXFEpQxxpQDDSvX44w2xwHw0pR3yMzO9Dmi/bMEZYwx5cTJrY8mrUoDNu5K5/3ZX/odzn5ZgjLGmHIiLiaWq7pfQEwghu8W/cKCTYv9DmmfLEEZY0w50qRqQ05pfRRBgrww5W2ycrL8DqlIlqCMMaacOb3NcdSvVId1Ozby0dxv/A6nSJagjDGmnImPjefq7hcQIMDXOpoVW1f7HVKhLEEZY0w51KJ6E45u3p+8YB4vT32fvGCe3yH9jSUoY4wpp85pfxJVKlRiYfpSflo6we9w/sYSlDHGlFPJCUlc2OlMAN6d9RnbMrf7HNHeLEEZY0w51qthFzrUac2urAzejrBukCxBGWNMORYIBLik8znEx8Qxdvkk5mxQv0P6kyUoY4wp5+qk1uLUNscCbnDD7NxsnyNyLEEZY4zh5FZHUi+1Nmt3bOArHe13OEAEJigRqSIiK0XkQu91gogMF5HNIrJJRP7rc4jGGFPmxMfGc2mXcwH4dN5I1u/c5HNEEZiggBeB+iGv7wUEaAZ0A4aIyAV+BGaMMWVZu9pCv8Y9yM7N5rVpHxAMBn2NJ6ISlIgMASoBs0MmDwHuV9UtqroceAy4wofwjDGmzBvc8TQqxicxc/08Jq6e7mssEZOgRKQJcDdwcci0KkBdYF7IoguA9qUbnTHGlA+VK1RiUIdTAXhj+sdkZPs3RHxEJCgRiQXeAW5W1fUhs1K83xkh0zKA5NKKzRhjypvDm/ahRfUmbMncxgc+jhsVEQkKuBNQVS34lNgu73dSyLRkYGepRGWMMeVQTCCGy7qcR0wghlGLx7Bk8wp/4vBlr393DnCGiGwVka24KrzngfuB9bhGEvlasXeVnzHGmGKWVrUBx7U8nGAwyMtT3yMvr/Q7k40r9T0WQlVbhb4WkZnAU6r6hojsBO4WkVm4Kr+bgad9CNMYY8qVs9oez28rp7F0y0q+XzKWY1oMKNX9R0oJal/uAuYAc4EpwKe4pujGGGNKUIX4ClzU+SwA3p/1BZt3by3V/UdECaogVe0Y8ncmMNT7McYYU4q61e9Al3rtmbZ2Nm/N+IR/9b601PYdDSUoY4wxPgkEAlzc+WwSYxOYsGoaM9eVXhMAS1DGGGP2qWbF6pzZ7ngAXp3+AVk5WaWyX0tQxhhj9uu4lgNpWLkeG3Zu4rP5o0pln5agjDHG7FdcTCyXdTkPgM8XjGLN9vX7WeOfswRljDEmLK1qNuPwpn3IzcvllWnvl3hnspagjDHGhG3QIaeQmpjC3I0LGbdiconuyxKUMcaYsKUmpjC4w2kAvDXzE3Zm7drPGgfPEpQxxpgD0j+tJ21qtmD7np28N+uLEtuPJShjjDEHJBAIcGnXc4mNiWX0knFszdxeIvuJyJ4kjDHGRLYGlepyVbfBzN+0mJSEiiWyD0tQxhhjDkq/tB70S+tRYtu3Kj5jjDERyRKUMcaYiGQJyhhjTESyBGWMMSYiWYIyxhgTkSxBGWOMiUjloZl5LMD69SXf864xxpgDE3Juji04rzwkqLoAgwYN8jsOY4wxRasLLAmdUB4S1BTgUGAdkOtzLMYYY/YWi0tOUwrOCJT0eB7GGGPMwbBGEsYYYyKSJShjjDERyRKUMcaYiGQJyhhjTESyBGWMMSYiWYIyxhgTkSxBGWOMiUiWoIwxxkSk8tCTRNhE5EjgIaAFsBF4VFVfEpEE4FngDFxvFE+o6oP+RXrgROQE4AGgCe7YHikrxwYgIlWAWcBdqvpGWTguEbkYeAnYEzJ5KPA+0X9sdYEXgMOATGC4qt4Z7Z+biAzCfWahkoAfgROI7mPrCTwDCLAJeEhVXynJz8wSlEdEGgKfAkOAL4AuwCgRWQ4MwH0ozYDKwHciskZV3/Il2APknQw+AU5V1ZEi0hn4VUSmAGcSxccW4kWgfsjre4n+4+oMPK6qt4ZOFJEHif5j+wKYBtTGdXMzRkTmA+2J4mNT1XeBd/Nfi0gn4HvgFqL4OykiMbjP7CZVfUdEugHjvHPIOZTQcVkV31/SgPdU9TNVzVPVKcAvQB9c0rpfVbeo6nLgMeAKvwI9UKq6DqjpJacYoDqQA+wgyo8NQESGAJWA2SGTo/64cBdJMwuZHtXHJiI9gKbAdaqaqarLcBeBPxPlxxZKROJxyeoeVf2d6D62qkAtICAiASCIO4dkUYLHZSUoj6qOA8blvxaRarhOZt/GXeHNC1l8Ae5KL2qo6g4RSQa24T73h3HF9Kg+NhFpAtwN9Aa+86ZVIfqPKxY4BBgsIk8AGcAruOqjqD42XOKdDdwjIhfiqvieB14l+o8t1FBgN/B8tH8nVTVdRJ4F3gRex3Xwej2uE+4SOy5LUIUQkcrAl8AkXDUEuBMEIX8nl3ZcxSATqIg78X2L++eBKD027yT+DnCzqq4XkfxZKd7vqDwuT01gKu6EcBrQGlfFkuDNj+Zjy7/4G4MrSbXCXVxs8uZH87EB4N2XuQW4UlWDIhLV30mv5iUTOA93K6Q3MALY6i1SIsdlCaoAEWmJOxHMAwbhbnAS8hvcm7+zlEP7x1Q1D1cknyoiw4Gu3qxoPbY7AVXVEQWm7/J+R+txoarrgf4hk2aKyP+AY73XUXtsuEYf21X1Hu/17yLyCq6qCKL72PIdA+QB33ivo/07eRrQR1Vv8V6PEZFXKeHPzO5BhRCRfrhS0+fAGV79+BZgPe4mYL5W7F2kjWgi0l9EphWYnAhE+7GdA5whIltFZCuuWuF54H6i+7gQkbYicm+ByQm4q9ioPjZcFVCyV8rIF0f0fx9DnQx85F0UUgbOIw1x54xQObhSb4kdl5WgPCLSDPgauF1V/1dg9tvA3SIyC1d9dDPwdCmH+E/MBOqLyI24uHsAlwCn4r5cUXlsqtoq9LWIzASe8pqZ7yRKj8uzFbhJRFbj7s10Aq4DrgHmEt3H9gPuxPa4iNyEO7ldAlwFLCW6jy1fT1wJP1Q0n0e+Bx4UkcuBl3EtTC8DLgVWUkLHZSWovwwFUnEfws6Qn4eBu4A5uBPDFFwd7Iv+hXpgVHUbcByumL4ZGA5cqqpjiPJj24eoPi5VXQOchGsNtR0X/zBV/YToP7ZMXPVlU9xN9u9wz+V9SpQfW4g0YG2BaVF7bKo6F3f+uAJ38fQecKuqfkEJHpeNqGuMMSYiWQnKGGNMRLIEZYwxJiJZgjLGGBORLEEZY4yJSJagjDHGRCRLUMYYYyKSJSgTkUQk6I1hFXEiLbbiiEdEnhaRi4qY187bR9o/2UdJEZGKIjJbRGr4HYspXpagjCnnRKQrMBDXMW3UUdVduF7eH/U7FlO8LEEZY+4GXszvNy5KvQqc7nVZZsoI64vPRCURORo3QmkH3OBpE4GrVXWBiAzADX6Xqqo7veXvAU5Q1a5eVdUy3GjCw3AdYU7BDY2g3vIdgMdxfaptxZ3A7wsJoauI3IUbumQxMNTrOqqwWFt42+qH6/V5EXCbqn7pzV8OPInrG7G7t73bVfUrb35VXNcxx+LG87oTNzZUc2+AuIL7uxm4Fjcw5UzccCQTi4itibfdK0Om1cSVSI7CdUX0ZIF1Ur3jOQP33v8EXK+qa8OJV0SCwH3A5bi+IDvjRmN9xnuPNgEfAneq6h5vmy33NV9Vd4vI97i+Cm8o7FhN9LESlIk6ItIYN17Xx0Ab4HDcGEMHWsVzD65vscNwg6495m2/Bu6kuxaXMC7Fddx6aci6V3nrt8d1cPqeN9JowVgDwFe40Yt7Ah1xg/W9XqA37//DndS74hJY6Pz3cf3WHQacD9yOGzDub0TkCtxAclfjOpj9FvjJS0SFOQ6Y5/X9l+9joA7Q1zvO/xRYZziug9ejcX3qBYFRIpJ/wRtOvINwn9sQIB4YhXsfOwGDccNVPO0dU4V9zQ/xHX8NR2LKACtBmWgUhysV5Pc6v0xE3sadlA/Efao6FkBEngdu9aafDWTjOtTNAuaJyNVAbsi6D6vqt966j+BGY64JbCywjyRc6eE1Vd3sLf8YbqiQ2sAqb7n3VfUDb/69wO9AmjcI49FAB1Wd5c2/FhhZxDHdhuvEM38coge8EuVQXC/TBXXFdfKJt+3WuKTTzusgFBH5Dy7pICJNvdgb5Cc1ERkM/AEcIyILw4x3eMj2L8K930NVNQioiFwJjBORfwOn72u+qm73tjnPbU5S8kvOJrpZgjJRR1WXiMgI78TZFjf+TAdcddSBWBjy93bclTy4UtksLznl7/P9AusuCfk7f1TRpALLoKoZIvICcJ7XGKElrkoL9i5VFIwFL57WuAH+ZofM/62wg/FGbW0EvCwiL4XMSvS2UZjaBY6lHbAnP3l4Jof83Sb/0EJGMAY3SJ0AFcKMN3SfbXElrh0h2wzganhahDE/f6yzdO93LaJnIECzD5agTNQRkXa4k94PuGHDX8GNcTXUW6SwLvoL+65nFXgdCJm+v27+cwuZVlgVX0XcIJi7gc9wVZM7gV/2E0v+9rIL224R8hPeEGBGgXm7i1gnr7Dti0jAK60UjC3Oi6kTf3+PNuOGcg8n3tB44nCfZ2HN3NeEMT9f/vEX9tmYKGT3oEw0GgLMUNXTVPVpr5quCXsnGIDKIes0PYDtLwTai0h+iQoRuVNEPjmIWAfgrvIPVdUHvKq32t68cE7kc3Ej6bYPmdatsAW9cb/WA/VVdXH+D67hwNFFbH89rmoy3yxciatjyLTOIX/Px5XsKoZsfx3u/l/LA4m3wDZbAKtDtlkNeMjb1v7m58t/Dmr9fvZnooSVoEwk6ywiOQWmzcZV5bQSkUNxV9An40b3zK/imYu7Qh8mIsNwN+uPxzU+CMe7uBaCz3v3i5riWobdeBDHkI47iZ4tIr/gTvZPefMKDqH9N6q6WES+wlXbDfW2lX/vrbBS3iPAXSKyDtcy8TxcgupfxC6m4RpV5O9PReRb4FWvwUWit83Q+V8Cb3nxbALuxzUAWaCqWw8wXoB3cIPevel9Xqm4UvFKVd0mIvucH7KdDsDs/JZ9JvpZCcpEsntxN9dDf47ENTcejWsdNw04AddMupaINPBuml+MK73MA07xthUWb/1jcfe2fse1rrtfVd840APwmnffjrvan4c70d4MbAG6hLmZi3GJeCzwEX89UFtYteDTuNaIj3j7Ows4Q1V/LWLb3wItRaROyLRzvXV/wjWOeKbAOkOAqcDnuCRYGThSVfPvxR1IvPkP2h4FVMXd7/qKv5LrfueH6Ad8XcRxmihkI+oaE8FEJBl3ch4Z8kxQN2A8rpqtYAnzYPbxDfCjqj5RDNsq8XiL2G9lYDVwiKouK4l9mNJnCcqYCCYiMcAG4C3gOaAKropwtaoWLEEc7D564Eo5bVX1HzUwKI14i9jvDUBHVR1SUvswpc+q+IyJYF73Qyfh7vHMxj2wOh/3gHFx7WMSrsq00M5iD3BbJR5vQV5LycuAW0pqH8YfVoIyxhgTkawEZYwxJiJZgjLGGBORLEEZY4yJSJagjDHGRCRLUMYYYyLS/wMYOBpvuj4I2QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot(sweep, color='C2')\n",
    "decorate(xlabel='Launch angle (degree)',\n",
    "         ylabel='Range (m)',\n",
    "         title='Range as a function of launch angle',\n",
    "         legend=False)\n",
    "\n",
    "savefig('figs/chap23-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `maximize` to search for the peak efficiently."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 degree 17.15514421874411 meter\n",
      "90 degree 1.2569667261582053e-14 meter\n",
      "34.376941012509455 degree 101.06319487191678 meter\n",
      "55.623058987490545 degree 93.24073676688748 meter\n",
      "21.246117974981075 degree 82.63917923456437 meter\n",
      "42.49223594996215 degree 103.32194321543766 meter\n",
      "47.507764050037856 degree 101.41737997289113 meter\n",
      "39.39246911258517 degree 103.27504916143039 meter\n",
      "44.40799721266088 degree 102.90451237123467 meter\n",
      "41.308230375283905 degree 103.41115569950708 meter\n",
      "40.576474687263406 degree 103.38512107429113 meter\n",
      "41.760480261941645 degree 103.4149581469886 meter\n",
      "42.03998606330441 degree 103.38949811498512 meter\n",
      "41.58773617664667 degree 103.41463212933446 meter\n",
      "41.867241978009446 degree 103.41261934653042 meter\n",
      "41.69449789271447 degree 103.41499899416574 meter\n",
      "41.65371854587384 degree 103.41492192196334 meter\n",
      "41.71970091510101 degree 103.41500757170726 meter\n",
      "41.73527723955511 degree 103.41499793513337 meter\n",
      "41.710074217168575 degree 103.41500782026304 meter\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x</th>\n",
       "      <td>41.710074217168575 degree</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fun</th>\n",
       "      <td>103.41500782026304 meter</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                         True\n",
       "x          41.710074217168575 degree\n",
       "fun         103.41500782026304 meter\n",
       "dtype: object"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bounds = [0, 90] * degree\n",
    "res = maximize(range_func, bounds, params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`res` is an `ModSimSeries` object with detailed results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x</th>\n",
       "      <td>41.710074217168575 degree</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fun</th>\n",
       "      <td>103.41500782026304 meter</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                         True\n",
       "x          41.710074217168575 degree\n",
       "fun         103.41500782026304 meter\n",
       "dtype: object"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`x` is the optimal angle and `fun` the optional range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "41.710074217168575 degree"
      ],
      "text/latex": [
       "$41.710074217168575\\ \\mathrm{degree}$"
      ],
      "text/plain": [
       "41.710074217168575 <Unit('degree')>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optimal_angle = res.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "103.41500782026304 meter"
      ],
      "text/latex": [
       "$103.41500782026304\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "103.41500782026304 <Unit('meter')>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max_x_dist = res.fun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Under the hood\n",
    "\n",
    "Read the source code for `maximize` and `minimize_scalar`, below.\n",
    "\n",
    "Add a print statement to `range_func` that prints `angle`.  Then run `maximize` again so you can see how many times it calls `range_func` and what the arguments are."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "def maximize_golden(max_func, bracket, *args, **options):\n",
      "    \"\"\"Find the maximum of a function by golden section search.\n",
      "\n",
      "    :param min_func: function to be maximized\n",
      "    :param bracket: interval containing a maximum\n",
      "    :param args: arguments passes to min_func\n",
      "    :param options: rtol and maxiter\n",
      "\n",
      "    :return: ModSimSeries\n",
      "    \"\"\"\n",
      "    def min_func(*args):\n",
      "        return -max_func(*args)\n",
      "\n",
      "    res = minimize_golden(min_func, bracket, *args, **options)\n",
      "\n",
      "    # we have to negate the function value before returning res\n",
      "    res.fun = -res.fun\n",
      "    return res\n",
      "\n"
     ]
    }
   ],
   "source": [
    "source_code(maximize)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "def minimize_scalar(min_func, bounds, *args, **options):\n",
      "    \"\"\"Finds the input value that minimizes `min_func`.\n",
      "\n",
      "    Wrapper for\n",
      "    https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize_scalar.html\n",
      "\n",
      "    min_func: computes the function to be minimized\n",
      "    bounds: sequence of two values, lower and upper bounds of the range to be searched\n",
      "    args: any additional positional arguments are passed to min_func\n",
      "    options: any keyword arguments are passed as options to minimize_scalar\n",
      "\n",
      "    returns: ModSimSeries object\n",
      "    \"\"\"\n",
      "    try:\n",
      "        min_func(bounds[0], *args)\n",
      "    except Exception as e:\n",
      "        msg = \"\"\"Before running scipy.integrate.minimize_scalar, I tried\n",
      "                 running the slope function you provided with the\n",
      "                 initial conditions in system and t=0, and I got\n",
      "                 the following error:\"\"\"\n",
      "        logger.error(msg)\n",
      "        raise(e)\n",
      "\n",
      "    underride(options, xatol=1e-3)\n",
      "\n",
      "    res = scipy.optimize.minimize_scalar(min_func,\n",
      "                                         bracket=bounds,\n",
      "                                         bounds=bounds,\n",
      "                                         args=args,\n",
      "                                         method='bounded',\n",
      "                                         options=options)\n",
      "\n",
      "    if not res.success:\n",
      "        msg = \"\"\"scipy.optimize.minimize_scalar did not succeed.\n",
      "                 The message it returned is %s\"\"\" % res.message\n",
      "        raise Exception(msg)\n",
      "\n",
      "    return ModSimSeries(res)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "source_code(minimize_scalar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Manny Ramirez problem\n",
    "\n",
    "Finally, let's solve the Manny Ramirez problem:\n",
    "\n",
    "*What is the minimum effort required to hit a home run in Fenway Park?*\n",
    "\n",
    "Fenway Park is a baseball stadium in Boston, Massachusetts.  One of its most famous features is the \"Green Monster\", which is a wall in left field that is unusually close to home plate, only 310 feet along the left field line.  To compensate for the short distance, the wall is unusually high, at 37 feet.\n",
    "\n",
    "Although the problem asks for a minimum, it is not an optimization problem.  Rather, we want to solve for the initial velocity that just barely gets the ball to the top of the wall, given that it is launched at the optimal angle.\n",
    "\n",
    "And we have to be careful about what we mean by \"optimal\".  For this problem, we don't want the longest range, we want the maximum height at the point where it reaches the wall.\n",
    "\n",
    "If you are ready to solve the problem on your own, go ahead.  Otherwise I will walk you through the process with an outline and some starter code.\n",
    "\n",
    "As a first step, write a function called `height_func` that takes a launch angle and a params as parameters, simulates the flights of a baseball, and returns the height of the baseball when it reaches a point 94.5 meters (310 feet) from home plate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def event_func(state, t, system):\n",
    "    \"\"\"Stop when the y coordinate is 0.\n",
    "    \n",
    "    state: State object\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: y coordinate\n",
    "    \"\"\"\n",
    "    R, V = state\n",
    "    return R.x - system.x_wall"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Always test the slope function with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "-94.5 meter"
      ],
      "text/latex": [
       "$-94.5\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "-94.5 <Unit('meter')>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "system = make_system(params)\n",
    "system.set(x_wall = 94.5 * m)\n",
    "event_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def height_func(angle, params):\n",
    "    \"\"\"Computes the height of the ball at the wall.\n",
    "    \n",
    "    angle: launch angle in degrees\n",
    "    params: Params object\n",
    "    \n",
    "    returns: height in meters\n",
    "    \"\"\"\n",
    "    params = Params(params, angle=angle)\n",
    "    system = make_system(params)\n",
    "    system.set(x_wall = 94.5 * m)\n",
    "\n",
    "    results, details = run_ode_solver(system, slope_func, events=event_func)\n",
    "    height = get_last_value(results.R).y\n",
    "    \n",
    "    return height"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your function with a launch angle of 45 degrees:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "10.953200948514532 meter"
      ],
      "text/latex": [
       "$10.953200948514532\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "10.953200948514532 <Unit('meter')>"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "height_func(45 * degree, params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use `maximize` to find the optimal angle.  Is it higher or lower than the angle that maximizes range?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>x</th>\n",
       "      <td>44.25082945182002 degree</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>fun</th>\n",
       "      <td>10.967930469934105 meter</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                        True\n",
       "x          44.25082945182002 degree\n",
       "fun        10.967930469934105 meter\n",
       "dtype: object"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "bounds = [0, 90] * degree\n",
    "res = maximize(height_func, bounds, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "44.25082945182002 degree"
      ],
      "text/latex": [
       "$44.25082945182002\\ \\mathrm{degree}$"
      ],
      "text/plain": [
       "44.25082945182002 <Unit('degree')>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "optimal_angle = res.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "10.967930469934105 meter"
      ],
      "text/latex": [
       "$10.967930469934105\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "10.967930469934105 <Unit('meter')>"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "optimal_height = res.fun"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With initial velocity 40 m/s and an optimal launch angle, the ball clears the Green Monster with a little room to spare.\n",
    "\n",
    "Which means we can get over the wall with a lower initial velocity."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Finding the minimum velocity\n",
    "\n",
    "Even though we are finding the \"minimum\" velocity, we are not really solving a minimization problem.  Rather, we want to find the velocity that makes the height at the wall exactly 11 m, given given that it's launched at the optimal angle.  And that's a job for `root_bisect`.\n",
    "\n",
    "Write an error function that takes a velocity and a `Params` object as parameters.  It should use `maximize` to find the highest possible height of the ball at the wall, for the given velocity.  Then it should return the difference between that optimal height and 11 meters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "def error_func(velocity, params):\n",
    "    \"\"\"Returns the optimal height at the wall minus the target height.\n",
    "    \n",
    "    velocity: initial velocity in m/s\n",
    "    params: Params object\n",
    "    \n",
    "    returns: height difference in meters\n",
    "    \"\"\"\n",
    "    print(velocity)\n",
    "    params = Params(params, velocity=velocity)\n",
    "    bounds = [0, 90] * degree\n",
    "    res = maximize(height_func, bounds, params)\n",
    "    return res.fun - 11*m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Test your error function before you call `root_bisect`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "40.0 meter / second\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "-0.03206953006589508 meter"
      ],
      "text/latex": [
       "$-0.03206953006589508\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "-0.03206953006589508 <Unit('meter')>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "error_func(40*m/s, params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then use `root_bisect` to find the answer to the problem, the minimum velocity that gets the ball out of the park."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30.0 meter / second\n",
      "50.0 meter / second\n",
      "40.0 meter / second\n",
      "45.0 meter / second\n",
      "42.5 meter / second\n",
      "41.25 meter / second\n",
      "40.625 meter / second\n",
      "40.3125 meter / second\n",
      "40.15625 meter / second\n",
      "40.078125 meter / second\n",
      "40.0390625 meter / second\n",
      "40.01953125 meter / second\n",
      "40.009765625 meter / second\n",
      "40.0048828125 meter / second\n",
      "40.00732421875 meter / second\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>converged</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>root</th>\n",
       "      <td>40.008544921875 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "converged                              True\n",
       "root         40.008544921875 meter / second\n",
       "dtype: object"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "bounds = [30, 50] * m/s\n",
    "res = root_bisect(error_func, bounds, params, rtol=1e-4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "40.008544921875 meter/second"
      ],
      "text/latex": [
       "$40.008544921875\\ \\frac{\\mathrm{meter}}{\\mathrm{second}}$"
      ],
      "text/plain": [
       "40.008544921875 <Unit('meter / second')>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "min_velocity = res.root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And just to check, run `error_func` with the value you found."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "40.008544921875 meter / second\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "0.002204447936565046 meter"
      ],
      "text/latex": [
       "$0.002204447936565046\\ \\mathrm{meter}$"
      ],
      "text/plain": [
       "0.002204447936565046 <Unit('meter')>"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "error_func(min_velocity, params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
